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Abstract. We present a numerical algorithm for the incorporation of the active cosmic ray transport, into the 
ZEUS-3D magnetohydrodynamical code. The cosmic ray transport is described by the diffusion-advection equa- 
tion. The applied form of the diffusion tensor allows for anisotropic diffusion of cosmic rays along and across the 
magnetic field direction, which is controlled by two parameters: the parallel and perpendicular diffusion coeffi- 
cients. The implemented numerical algorithm is tested by comparison of the diffusive transport of cosmic rays 
to analytical solutions of the diffusion equation. Our method is numerically stable for a wide range of diffu- 
sion coefficients, including the realistic values inferred from the observational data for the Milky Way of about 
6 X 10^* cm^ s~^. The presented algorithm is applied for exemplary simulations of the the Parker instability 
triggered by cosmic rays injected by a single SN remnant. 
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1. Introduction 

One of the major components of the interstellar medium 
(ISM), the cosmic ray (CR) gas, consists of relativis- 
tic electrons, protons and heavier atomic nuclei (see eg. 
Berezinski et al. 1990). It was shown beyond any doubt 
that the cosmic ray particles are accelerated in the process 
of diffusive acceleration by shocks associated with super- 
nova remnants (SNR) in galactic disks (e.g. Koyama et al. 
1995). Recent models suggest that the conversion rate of 
the supernova energy into cosmic ray energy is in the range 
of 10 - 50 % (see eg. Jones 1998 and references therein). 
The total kinetic energy output from a single supernova 
is of the order of 10^^ erg, therefore the total CR energy 
per unit volume, produced within a supernova remnant, is 
significant as compared to thermal, kinetic and magnetic 
energy densities of the ISM. 

Although the velocity of individual CR particles is 
close to the speed of light, the bulk motion of CR is dif- 
fusive and the CR bulk speed is of the order of Alfven 
speed, i.e. typically a few tens of km/s. Recent studies 
by Giacalone and Jokipii (1999) and Jokipii (1999) sug- 
gest that the diffusion of cosmic ray gas in a turbulent 
magnetic field proceeds preferentially along the direction 
of the mean magnetic field. In our case the term cosmic 
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rays means protons and nuclei but not electrons, since 
their contribution to the pressure is negligible. The esti- 
mations made by Strong and Moskalenko ( 1998) based on 
the GALPROP model provide the parallel diffusion coef- 
ficients of the order of = 6 x 10^* cm^ s~^. This value 
is 2-3 orders of magnitude larger than the diffusion coeffi- 
cient for turbulent mixing of the ISM. The large energies 
carried by the cosmic ray component as well as its highly 
diffusive nature imply that the cosmic ray component can- 
not be neglected in the studies of dynamics of the ISM. 
That statement follows directly from investigations of sta- 
bility of the ISM on spatial scales of the order of one up 
to a few kiloparsecs. Parker (1966, 1967) found that the 
multicomponent interstellar medium stratified by vertical 
gravity is subject to an instability which is caused by the 
buoyancy of the weightless ISM components, i.e magnetic 
field and cosmic rays. 

The Parker instability has been extensively studied by 
numerous authors in the linear approximation under var- 
ious circumstances like different disk gravity models (Giz 
& Shu 1993; Kim & Hong 1998), rigid and differential ro- 
tation (Shu 1974; Fogfizzo & Tagger 1994, 1995; Hanasz & 
Lesch 1997) the presence of random magnetic field com- 
ponent (Parker & Jokipii 2000; Kim & Ryu 2001) and 
nonadiabatic effects in the ISM (Kosinski & Hanasz 2003). 
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Majority of the work was done within the hmit of very 
large diffusion of cosmic rays along magnetic field lines and 
negligible diffusion across magnetic field lines. The effect 
of finite diffusion was studied by Kuznetsov and Ptuskin 
(1983) and recently by Ryu et al. (2003), who demon- 
strate that within the linear approximation incorporation 
of the diffusion-advection equation and realistic diffusion 
coefficients leads to results consistent with the mentioned 
simplified description of cosmic ray transport. According 
to the analysis done by Ryu et al. (2003) the finiteness of 
the diffusion coefficient decreases the growth rate of the 
Parker instability. 

On the other hand numerical studies of the Parker in- 
stability investigate the effects of uniform vertical gravity 
(Kim et al. 1998), realistic vertical gravity (Kim et al. 
2000), selfgravity (Chou et al. 2000), the effects of spiral 
arms (Franco et al. 2002), finite resistivity (Hanasz et al. 
2002, Tanuma 2003; Kowal et al 2003), partial ionization 
(Birk 2002) and coupling to other disk instabilities (Kim, 
Ostriker & Stone 2002). 

Surprisingly, the powerful cosmic ray component, 
which according to the linear analysis is crucial for the 
growth rate of the Parker instability, is neglected in nu- 
merical studies except the recent paper by Hanasz & Lesch 
(2000) who incorporate the diffusion-advection equation 
for the cosmic ray transport along magnetic field and 
study the Parker instability triggered by cosmic ray in- 
jection in SN remnants, applying the thin fluxtube ap- 
proximation. That paper demonstrates the importance of 
the cosmic ray component for the global dynamics of the 
ISM, including the hydromagnctic dynamo effect. 

In this paper we describe how to introduce 
the cosmic ray component, within the diffusion- 
advection equation, into the ZEUS-3D MHD code 
(Stone & Norman 1992a, 1992b) developped at the 
Laboratory of Computational Astrophysics (NCSA, 
University of Illinois at Urbana Champaign, see 
http:/ /lea, ncsa. uiuc. edu/lca- codes-docs.htm^ . The ZEUS 
code uses a time-explicit, operator-split, finite-difference 
method to solve the MHD equations on a staggered mesh. 
The MHD algorithm employs the constrained transport 
formalism and the method of characteristics for accu- 
rate propagation of Alfven waves (Evans & Hawley 1988, 
Hawley & Stone 1995). 

In the present paper we focus on the numerical method 
for the active cosmic-ray transport. In Section 2. we intro- 
duce the set of basic equations. The numerical algorithm 
is described in section 3, followed by tests of the numerical 
method and a comparison of some results of computations 
to analytical solutions in section 4. Section 5 contains as 
an example an application of the extended code for stud- 
ies of the Parker instability triggered by the cosmic ray 
injection in a single SN remnant. Finally in Section 6. we 
summarize our results. 



2. Equations of MHD including cosmic-ray 
transport 

The diffusive cosmic ray (CR) transport on macroscopic 
astrophysical scales is described by the diffusion-advection 
equation. Following Schlickeiser and Lerche (1985) we ap- 
ply the following form of the transport equation 

dec 



dt 



V {e^rv) = V KVCcr ]-Pcr{V ■v) + Q 



(1) 



where Ccr and Pcr — (7cr — l)ecr are the cosmic ray energy 
density and cosmic ray pressure, 7cr (=4/3 in this paper) 
denotes the ratio of the specific heats of the relativistic 
cosmic ray gas, K presents the diffusion tensor, v is ve- 
locity of the thermal gas and Q is the source term for the 
CR energy density resulting from the cosmic ray injection 
by supernova remnants (SNR) or alternative sources. We 
note that eq. Q assumes that cosmic rays are treated as a 
magnetized relativistic gas. This assumption holds as long 
as the particles are tied to the magnetic field, i.e. as long 
as they gyroradius is significantly smaller than the char- 
acteristic spatial scales of the magnetic field. Only cosmic 
rays with ultrahigh energies are ruled out by our approach 
since their gyroradius is larger than the thickness of the 
galactic disk. Eq. |^ implicitly assumes that the diffu- 
sion tensor describes the interaction of charged particles 
with magnetic fiuctuations which appear on spatial scales 
considerably smaller than the characteristic scales in the 
interstellar medium. In the limit of low energies (< 100 
MeV) the the current approximation is not valid because 
cosmic-ray energy losses become important, although, the 
cosmic ray pressure still resides in that range. 

In the present approach we apply the concept of 
anisotropic cosmic ray diffusion following Giaccalone and 
Jokiph (1999), Jokipu (1999), Hanasz and Lesch (2000) 
and Ryu et al (2003). In order to describe formally the 
anisotropic cosmic ray diffusion we implement the diffu- 
sion tensor (see e.g. Ryu et al. 2003) of the form 



Kij = K^Sij + (i^ii - K^_)ninj, 



(2) 



where rii = Bi/ B are components of the unit vectors tan- 
gent to magnetic field lines. The above cosmic-ray trans- 
port equation supplements the standard set of ideal MHD 
equations 



| + V.(p.)^0, 

— + V-[ev) = -p{V-v), 

— + {VV)V^ V{p+pcr + — 

ot p \ 8n 



B ■ VB 



dB 

'dt 



V X (t) X J5) 



(3) 
(4) 
(5) 
(6) 



where the gradient of cosmic ray pressure Vpcr has been 
included in the equation of gas motion (see Berezinski et 
al. 1990). The other symbols have their usual meaning. 
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3. Numerical algorithm for the cosmic ray 
transport 

The actual form of the diffusion-advection equation Q 
is similar to the energy equation except the diffusion 
term, therefore we incorporate an integration algorithm 
for the advection part of the cosmic-ray transport, follow- 
ing the method of integration of the energy equation (see 
Stone and Norman 1992a,b). 

The integration method for the energy equation con- 
sists of a source step and a transport step. In the source 
step the — pVv term is evaluated together with possible 
explicit sources of the internal energy. In the transport 
step fluxes of internal energy, through cell boundaries, are 
computed using directional splitting. The total amount of 
internal energy within each cell is subsequently updated 
according to the sum of fluxes through all cell boundaries. 

The implementation of the cosmic ray transport re- 
quires an additional contribution of diffusive fluxes 



(7) 



corresponding to the term V yK^Ccrj in the cosmic ray 
diffusion-advection equation. The tcnsorial form of the dif- 
fusion coefficient K is used to describe the anisotropic cos- 
mic ray diffusion. 

In order to incorporate the diffusion of cosmic rays 
in the numerical algorithm, along the magnetic field lines 
one should compute first components of the unit vector 
n = B/B parallel to the magnetic field direction, sep- 
arately for each cell face. Since in the ZEUS code vector 
field components are centered on different cell faces an av- 
eraging is necessary for these vector components which are 
parallel to the given cell face. For instance, the magnetic 
field components on 1-faces (assigned with the superscript 
'If') are given by 

^llj,k) = [Bi(ij,k}, (8) 

0.25(i32(ij,fc) + -B2(i-lj,fe) + 52(!-l,i+l,fc) + -S2(ij + i,fe)), 
0.25(-B3(ij^fc) + i?3(i_ij,fc) + ^3(i-lj,fe+l) + -S3(ij^fe+i))]. 

The field components on the other faces are computed 
analogously. 

The next step is a computation of cosmic-ray diffusive 
fluxes across cell interfaces. This requires a prior compu- 
tation of components of the gradient of cosmic-ray en- 
ergy density. All three gradient components contributing 
to fluxes through a given cell-face should be centered at 
the center of that cell-face. Moreover, a monotinization 
of derivatives is essential for the numerical stability of the 
overall algorithm as soon as cosmic rays are coupled to the 
gas dynamics through the Vpcr term in the equation of 
gas motion. We apply the following formulae for numeri- 
cal derivatives needed for the the flux components through 
the 1-faces: 

(Vecr)(ij_fc) = [{dxecr){i,j,k), idyecr){i,j,k), idzecr){i,j,k)], (9) 

where 

{dxecr){i,j,k) = (ecr(ij,fe) - ecr{i-l.j,k))/ix(i) - (10) 



X(l-|- Sign(l,dy_iecr dy.rBcr)), 



(11) 



(rfzecr)(ij\fc) = 0.25(4,/ecr -I- dz^rCcr) (12) 
X(l -|-Sign(l,d^^;ecr dz^rCcr)), 

and the left and right derivatives used in the above for- 
mulae are given by 

dy^lBcr = 0.5((ecr(i-lj,fe) + ecr(ij,fc)) (13) 
-(ecr(i-lj-l,fc) + ecr(ij-l,fc)))/(2/(j) ~y{j-i)), 



— 0.5((ecr(i-lj + l,fe) + ecr(ij + l,fc)) 



(14) 



-(ecr(i-ij,fc) + ecr(i,i,fc)))/(?/(j+i) -yu)), 

dz^lCcr = 0.5((ecr(i-lj,fc) + ecr(ij,fc)) (15) 
— (ecr(i-lj-,fc-l) + ecr(i,j,fc-l)))/(^(/£) " ^(fc-l))i 

dz.rScr = 0.5((ecr(i_ij,fc+i) + eci.(ij,fc+l)) (16) 
-(ecr(i-lj,fc) + ecr(ij,fc)))/(^(fe+l) - Z{k))- 

The monotonized derivatives used in formulae (|ll|l and 
(fT^ reduce to centered derivatives if signs of left and right 
derivatives are the same. The cosmic ray energy fluxes 
through the cell-faces in the y and z- directions are con- 
structed in the analogous way. The cosmic ray diffusive 
fluxes through each cell-face can be now computed ac- 
cording to the formula Q. 

The standard stability analysis (see eg. Fletcher 1991) 
imposes the following necessary stability condition for the 
explicit numerical solutions of the diffusion equation: 

(min(Aa;, Ay, Az)-^ 



At < Co 



K 



(17) 



where Ccr — 0.5 is the Courant number corresponding to 
the diffusion problem. Numerical tests show that a slightly 
lower value, namely 0.3 suits better for the current prob- 
lem. 

The above limitation for the timestep, which is 
quadratic in the cell size, together with the application of 
monotonized derivatives implemented according formulae 
(|ll|l and H12() . ensures a stable numerical scheme for the 
active cosmic ray transport in the ZEUS code. 

The standard types of boundary conditions for the 
cosmic ray component can be implemented in the sim- 
ilar way as for the internal energy, except the outflow 
boundary condition. In case of the internal energy the 
outflow boundary condition relics on the replication of the 
contents of the starting and ending cells to the adjacent 
ghost zones. This kind of the outflow boundary condi- 
tion is, however, not appropriate for the diffusing cosmic 
ray component, since replication of cell contents means 
no gradients across the boundary. If gradients of cosmic 
ray energy vanish then only advection and no diffusion is 
possible across the boundary. In order to make it possi- 
ble for cosmic rays to diffuse across the boundary of the 
physical domain we first compute gradients of cosmic ray 
energy density across the boundary and then perform a 
linear extrapolation of the cosmic ray energy density from 
the interior cells to the ghost zones. 
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4. Test problems 

In this section we present several tests of our new nu- 
merical algorithm for the cosmic-ray difFusion-advection 
problem. 

In the following considerations we apply units which 
are convenient for the investigations of the dynamics of 
ISM on large spatial scales. The unit of length and time 
are 1 pc and 1 Myr, respectively. The unit of velocity is 

1 pc Myr~^ ~ 1 km s~^. The density is given in cm~^ 
corresponding to the number density of hydrogen atoms. 
The unit of the magnetic field is 1/iG. The diffusion co- 
efficient is expressed in units of 1 pc^ Myr^^. The realis- 
tic parallel diffusion coefficients for the cosmic ray trans- 
port as estimated by Strong and Moskalenko (1998), is 
based on the GALPROP model: K\\=Qx lO^^ cm^ s^^ = 

2 X 10^ pc^ 



Myr"\ In our numerical tests we apply diffu- 
sion coefficients ranging from 10^ 10^ pc^ Myr~^. 

We present test problems which are performed in the 
computational box of physical size 1000 pc x 1000 pc x 
500 pc with a spatial resolution of 100 x 100 x 50 grid 
zones in the x, y and z directions, respectively. Periodic 
boundary conditions are applied on all the domain bound- 
aries. 

4.1. Passive cosmic ray transport in one dimension 

As a first step we perform a test of passive cosmic ray 
transport along the magnetic field, which is directed along 
the X-axis, assuming K± — and a static medium, i.e 
V = 0. The passive transport means that the cosmic ray 
gas has no dynamical influence onto the motion of ther- 
mal gas. This can be achieved by neglecting the Vpcr m 
the equation of gas motion. In case of a static distribu- 
tion of gas and a uniform magnetic field parallel to the 
X-axis the cosmic ray transport equation reduces to a one- 
dimensional diffusion problem, described by the equation 



t = O.OOE + 00 



100 200 300 400 500 
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t=l .OOE + 02 
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Fig. 1. Comparison of numerical solutions of the ID, pure 
diffusion problem to analytical solutions for t = (upper 
panel) and t — 100 (lower panel). Continuous lines denote 
analytical solutions for different times and crosses repre- 
sent numerical solutions. 



dec 



dt 



K 



^j^g^ 4.2. Passive cosmic ray transport along an inclined 
magnetic field 



Assuming that the initial condition is given by the cosmic- 
ray distribution 



Ccro = A exp 



(19) 



where rg denotes the initial half-width of the Gaussian 
profile, we expect that the numerical solution should be 
close to the following analytical solution at any time 



,{x,t)=A^ 



AKt 



exp 



AKt 



(20) 



The comparison of the numerical solutions with the ana- 
lytical solution (eq. 20) is shown in Fig. ^ A perfect con- 
sistency of the both analytical and numerical solutions is 
obvious. 



The next test for the diffusive cosmic-ray propagation is 
to check if the simulated diffusion follows the analytical 
solution in case of an inclined magnetic field. We set up the 
values of if || — 100 and K± = and perform simulations 
for B^^By^O and B, = 0. 

Fig.|21presents two snapshots of cosmic ray distribution 
at t = and t = 100. It is apparent that qualitatively the 
anisotropic transport of cosmic ray energy proceeds along 
the magnetic field according to our expectations. 

However, a more precise evaluation of the numerical al- 
gorithm can be performed by fitting a 2D-Gaussian func- 
tion for the distribution of cosmic ray energy density in 
the xy-plane. The diffusion along the z-direction is absent 
in the current setup. The fitting procedure is performed 
with the aid of IDL routine gauss2dfit, which returns pa- 
rameters of a 2D-Gaussian distribution, i.e the semi-axis 
of the ellipsoid, as well as the amplitude of the peak and 
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Fig. 2. Diffusion of cosmic rays along an inclined magnetic 
field: (a) the initial spheroidal distribution of ecr at t = 
and (b) the ellipsoidal distribution at < = 100. 



the level of the constant background distribution of cosmic 
ray energy. The results of the fitting are shown in Fig. |31 
where two 1-d slices along the major and minor axis of the 
ellipsoid (crosses and asterisks respectively) are shown to- 
gether with the corresponding curves calculated on the 
base of the derived parameters of the ellipsoid. 

In the present case the fitted values of the parallel and 
perpendicular widths of the 2D Gaussian profile and the 
fitted amplitude are ry/ = 206.429, rj./ = 54.705, Af = 
11.096. Since K± = 0, the corresponding parameters of 
the exact solution can be derived from formula (|20l) repre- 
senting the time evolution of ID Gaussian profile. The ex- 
act analytical solution gives = \/rQ + AKt = 206.155 
and Aa = A/ yjr^ -f AKt = 11.884. If no perpendicular 
diffusion is present, the perpendicular widths should not 
change, i.e. r^a — 50.000. 

We note that the parallel width of the 2D Gaussian 
profile of the simulated cosmic ray distribution coincides 



t=1 .OOE + Oi 




00 



200 300 

X 



400 



500 



Fig. 3. Diffusion of cosmic rays along an inclined mag- 
netic field sX t — 100. Two cuts of the ellipsoid shown in 
Fig. [3 are made (b) along the major axis (crosses) and 
along the minor axis (asterisks). The full and dotted lines 
represent the corresponding cuts of the fitted 2D Gaussian 
profile. 



well with the analytical solution. However, there are two 
noticeable effects of the limited accuracy of our numerical 
algorithm. The first effect is the broadening of the per- 
pendicular widths r_L by about 10% of the original value. 
The broadening leads to a lowering of the peak amplitude 
Af with respect to the exact solution since the total cos- 
mic ray energy is conserved in absence of gas flows. The 
second effect of the limited numerical accuracy affects the 
formation of regions with negative values of cosmic ray 
energy density at the base of steep sides of the ellipsoid. 

The depth of the Ccr deficit (currently equal to 
—0.09123 aX t — 100) is strongly dependent on the grid 
resolution and the steepness of the initial cosmic ray en- 
ergy distribution across the magnetic field lines. The pres- 
ence of negative CR energy density regions may lead to 
significant numerical artifacts as soon as cosmic rays are 
coupled to the equation of gas motion. Therefore the spa- 
tial resolution of the grid in conjunction with the mag- 
nitude of the cosmic ray gradients is an issue of primary 
importance. The deficit vanishes in proportion to the grid 
resolution, however there is a limited freedom of reducing 
the size of the grid cell due to a drastic reduction of the 
timestep given by formula (|17|l . 

An alternative way is to incorporate an explicit per- 
pendicular diffusion given by equal to a few up to a few 
10 percent of K\\ . This procedure can be physically justi- 
fied by studies of cosmic ray transport in turbulent mag- 
netic fields (Giaccalone and Jokipii 1999). Another way to 
eliminate the spots of negative cosmic ray energy density 
is to add a smooth background of cosmic rays prior to the 
injection of very localized portions of cosmic ray energy. 
This is also a physically justified procedure, since e.g. the 
smooth background of cosmic rays is present in the ISM. 
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4.3. Active cosmic ray transport 

After testing the passive diffusion of cosmic rays we can 
now describe the active propagation of cosmic rays within 
the overall setup similar to that of the previous subsec- 
tion. In the present configuration we apply the uniform 
background of cosmic rays resulting from the assumption 
that cosmic ray pressure is equal to the magnetic and gas 
pressure 



t= 40 z= 



Ptot = Pgas + Pmag + Pc 

where 

Pmag ^pPgas 
Per — PpPgas 



(21) 

(22) 
(23) 



and ttp and /3p are in general constants of the order of 1 
in the ISM. In the present case we adopt ap ~ f3p ~ 1. 
The constant scaling factors between cosmic ray and gas 
background pressures as well as between magnetic and 
gas pressures is adopted for testing purposes only. This 
approach is useful especially in linear studies of Parker 
instability and follows the original work by Parker (1966, 
1967). In general, one can start with any arbitrarily chosen 
initial magnetic field and cosmic ray distributions. The 
corresponding background cosmic ray energy density in 
our units is Ccr — (7cr — l)^^apCgiPo — 147, where Cgi = 
7 km s""'^ denotes the isothermal sound speed and po = I 
is the background gas density. 

In order to investigate the coupling of cosmic rays 
to gas and magnetic field we switch on the Vpcr term 
in the equation of gas motion and inject about 50 % of 
20^1 erg of the kinetic output of the SN explosion. In our 
scaled units this corresponds to the initial peak amplitude 
of the Gaussian CR distribution equal to 100 times the 
background thermal energy density, i.e lOOCgjPo — 4900. 
Currently we apply a rather small diffusion coefficient of 
cosmic rays i^n — 100 (corresponding to 3 • 10^'"' cm^ s^^) 
and K± = 4 in order to illustrate the effects of cosmic ray 
propagation qualitatively. 

In Fig.^we present the cosmic ray - magnetic field dis- 
tribution (panel a) and the density - velocity distribution 
(panel b) at t = 40. The difference between the passive 
and active cosmic ray transport is remarkable. First of all 
the gradient of the cosmic ray pressure leads to the accel- 
eration of gas. Due to the effect of magnetic tension gas 
accelerates preferentially along the magnetic field up to a 
few km s~^. The outflow of gas from the cosmic ray in- 
jection region together with the expansion forced by the 
cosmic ray pressure implies the formation of a cavity in 
the gas distribution. On the other hand, gas accumulates 
outside the cavity as it is visible in^(b), in the upper-left 
and lower-right corners of the graphic. 

The gas motion along the magnetic field lines leads to 
an advection of cosmic rays with the gas flows, which is 
noticeable in Fig.^(a) as an enhanced cosmic ray energy 
density coinciding with the enhanced gas density in Fig. 0] 
(b). The coupling of cosmic rays to the gas component 
implies that the cosmic rays spread faster with respect to 
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Fig. 4. Diffusion of cosmic rays along an inclined magnetic 
field at i = 40: (a) the cosmic ray energy density ecr and 
magnetic field and (b) distribution of gas density and gas 
velocity. 

the passive (pure diffusion) transport. A broadening of the 
cosmic ray profile across the magnetic field is due to the 
pressure of the cosmic ray gas and partially due to the 
imposed perpendicular diffusion. 

5. Active cosmic ray transport in a vertically 
stratified atmosphere 

For the simulations of the cosmic ray transport in a strati- 
fied atmosphere we adopt a physical domain and grid sizes 
500 X 1000 X 1000 pc and 50 x 100 x 100 grid zones, in x, y 
and z directions respectively. We apply periodic boundary 
conditions to all the vertical domain boundaries, a reflec- 
tion boundary condition to the lower domain boundary 
and outflow condition to the upper boundary. 

The goal of the present work is to incorporate the cos- 
mic ray transport into studies of the dynamics of a grav- 
itationally stratified interstellar medium. In this section 
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Fig. 5. Geometry of the initial state stratified by vertical 
gravity: (a) the slice in yz-plane showing the stratification 
of the background distribution of cosmic rays and mag- 
netic field along with the localized cosmic ray injection 
region, (b) the initial, horizontal magnetic field is parallel 
to the horizontal and parallel to the y-axis. 

we perform an experiment similar to the ones presented 
in the previous sections. However, in the present case a 
uniform vertical gravity is taken into account for the con- 
struction of an initial equilibrium state. The equilibrium 
fulfills the magnetohydrostatic force balance equation 



dpt, 



dz 



0, 



(24) 



where ptot = pi^ + Up + (3p) denotes the total pressure 
and gz = —0.65 pc Myr~^ is the vertical, uniform gravi- 
tational acceleration. The center of cosmic ray injection is 
placed at X = 0, y = and z — 100. Two slices illustrating 
the geometry of the initial state are shown in Fig. |S1 

Fig. shows the state of the system at i = 100 in 
case of K\\ = 10'' (corresponding to 3 • 10^^ cm^ s~^) and 
K^^ — 0. Cosmic rays injected into a localized region dif- 
fuse anisotropically along the magnetic field lines and pop- 
ulate a fluxtube marked by magnetic lines threading the 
initial injection volume. Due to an excess of cosmic ray 
pressure the flux tube becomes underdense and its cen- 
tral part starts to rise against vertical gravity. The overall 
evolution of the fluxtube follows closely the one described 
in the thin fluxtube approximation by Hanasz & Lesch 
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Fig. 6. Propagation of cosmic rays in a vertically strat- 
ified atmosphere for the diffusion coefficient K\\ — 10^, 
after the local injection of cosmic ray energy: (a) cosmic 
ray energy density and magnetic field in the yz-plane at 
X = and (b) in the xy-plane at z = 400, (c) the density 
perturbation Ap/po and the gas velocity in the xz-plane. 



(2000). The gradient of the cosmic ray pressure acceler- 
ates gas, along the direction of magnetic field, reducing 
additionally the gas density at the neighborhood of the 
injection region. That effect enhances the strength of the 
buoyancy force. 
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Fig. 7. Same as Fig.lHlfor the diffusion coefficient = Fig. 8. Same as Fig.Qat t = 150. 
10^ at t = 100. 



At t — 100 the cosmic ray populated flux tube forms a 
rising Parker loop. In Fig.El(c) the relative density Ap/po 
is shown in the xz-plane together with the velocity field. 
The apparent tube-like cavity in the density distribution 
results from the local excess of cosmic rays. The upward 
gas velocity is a consequence of the buoyancy force. The 
rising tube compresses the overlying gas and pushes it 
toward higher altitudes. The system perturbed by cos- 
mic rays injected in a localized spherical region forms a 



buoyant fluxtube and evolves in a fashion resembling the 
development of an undulatory Parker instability mode. 

When the diffusion coefficient takes a realistic value 
i^il = 10^, which corresponds to 3 ■ 10^^ cm^ s~^ the evo- 
lution of the system is remarkably different (see Figs. 
and IHl The distribution of cosmic rays along the flux 
tube becomes relatively uniform before the buoyancy force 
starts to displace the tube in the vertical direction. The 
perturbation provided by the cosmic ray input excites ini- 
tially (up to < = 100) the interchange mode of the Parker 
instability with a weak contribution of the undulatory 
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mode. Later on, at t = 150, the growing contribution of 
the undulatory mode becomes apparent. 

Due to the more efficient diffusion, cosmic rays fill in 
initially a larger volume, compared to the case of lower 
diffusion coefficient. Therefore the excess of cosmic ray 
pressure, and hence the buoyancy force, becomes weaker 
but distributed over a larger volume. At the fixed time 
t — 100 after the cosmic ray injection the maximum verti- 
cal speed is smaller compared to the case of smaller diffu- 
sion coefficient i^n — 10^, although later on the instabil- 
ity accelerates following the emergence of the undulatory 
mode. The apparent tendency seems to be opposite to that 
resulting from the linear stability analysis by Ryu et al. 
(2003). However, we point out that the lower values of the 
diffusion coefficient are clearly leading to stronger nonlin- 
ear effects, (for instance the vertical speed of the buoy- 
ant gas is almost equal to the sound speed). The large 
cosmic ray energy remains localized in a more compact 
region, therefore the applicability of the linear approxi- 
mation for the discussed case of lower parallel diffusion 
coefficient (i^n = 10^) is questionable. Despite the men- 
tioned differences the late stages of the system for small 
and realistic diffusion coefficients remains rather similar. 



6. Summary and conclusions 

In this paper we presented a numerical algorithm for 
the inclusion of cosmic ray dynamics, described by the 
diffusion- ad vection equation, into the MHD code ZEUS- 
3D. In order to check the presented method we compared 
results of the diffusive passive transport of cosmic rays 
with analytical solutions of the diffusion equation. Our 
method appeared to be numerically stable in case of ac- 
tive transport for a wide range of diffusion coefficients, 
including the realistic values inferred from the observa- 
tional data by Strong & Moskalenko (1998) for the Milky 
Way. 

We applied the presented numerical algorithm to two 
exemplary simulations of the excitation of the Parker in- 
stability triggered by cosmic rays injected by a single SN 
remnant. The only difference between the input parame- 
ters of the two simulations is the magnitude of the par- 
allel diffusion coefficient. The simulation corresponding 
to the realistic value of the parallel diffusion coefficient 
if II = 3 X 10^^ cm^ s^^, presented in Figs. [3 and |H1 ap- 
peared to develop the Parker instability mode slower than 
the one performed for K^^ — 3 x 10^^ cm^ s~^, presented 
in Fig. El Such a tendency differs from the one resulting 
from linear analysis of the Parker instability by Ryu et 
al. (2003). We show these two examples of evolution of 
the system in order to demonstrate that in some circum- 
stances the finiteness of the diffusion coefficient may lead 
to effects which can not be described within the simple lin- 
ear approximation. Therefore a verification of all former 
analytical and numerical results concerning the nonlinear 
development of the Parker instability in presence of cosmic 
rays is necessary. 



The presented work is just a starting point, which fo- 
cuses on developing the basic computational techniques. 
In the next step we plan to combine the cosmic ray trans- 
port in a more realistic application by including the dy- 
namo action of the cosmic ray component, reconfiguration 
of the magnetic field by magnetic reconnection, different 
random spatial cosmic ray source distributions and differ- 
ent supernova rates. The future work should also include 
effects of cosmic ray losses and extensions of the present 
algorithm to the energy dependent cosmic ray transport. 
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